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Executive  Summary 

Maximum-likelihood  (ML)  techniques  are  useful  in  finding  synchronizer  structures  for  various  cases.  Synchronizers 
for  frequency,  phase,  and  timing  have  been  found  for  various  bandpass  signaling  techniques  such  as  PSK,  DPSK,  QAM, 
MSK,  and  CPM.  These  include  data-aided,  decision-directed,  and  clock-aided  cases.  This  report  describes  a  new  non- 
data-aided,  non-decision-directed  ML-based  phase  and  timing  synchronizer  for  a  dual-h,  full-response,  4-ary  CPM 
waveform.  The  derived  structure  will  be  incorporated  into  future  simulations  to  compare  performance  among  several 
possible  phase  and  timing  synchronizers. 
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I.  Introduction 

Successful  operation  of  digital  communication  systems  requires  that  receivers  achieve  synchronization.  Carrier 
frequency  and  phase,  waveform  time-of-arrival,  frame  synchronization  for  time-division  multiplexing  (TDM)/time- 
division  multiple  access  (TDMA),  and  timing  for  a  frequency-hopping  pattern  or  for  a  direct  sequence  spreading  code 
for  spread-spectrum  are  examples  of  typical  unknowns  that  must  be  estimated  for  good  performance. 

There  are  many  successful  synchronization  schemes.  Mengali  and  D’ Andrea  [1]  present  theoretical  and  practical  de¬ 
tails  for  synchronizers  of  frequency,  phase,  and  symbol  timing  for  various  baseband  and  bandpass  signaling  techniques, 
including  phase-shift  keying  (PSK),  differential  phase-shift  keying  (DPSK),  quadrature  amplitude  modulation  (QAM), 
minimum-shift  keying  (MSK),  and  continuous-phase  modulation  (CPM).  This  paper  generalizes  their  joint  phase  and 
timing  synchronizer  results  for  CPM  with  a  single  modulation  index  (single-/)  CPM)  to  CPM  with  two  modulation 
indices  (dual-/)  CPM),  assuming  a  full-response,  rectangular,  instantaneous  frequency  pulse  (1REC)  and  a  4-ary  sig¬ 
naling  alphabet  [2].  In  addition  to  frequency,  phase,  and  symbol  timing,  the  dual-/)  case  also  requires  the  so-called 
super-baud  synchronization  inherent  with  CPM  using  more  than  one  modulation  index.  A  previous  report  by  Pettit 
and  Wahlen  [3]  presented  a  frequency  synchronizer  for  the  same  case. 

Synchronizers  are  developed  using  either  an  ad  hoc  approach,  based  upon  intuitive  reasoning,  or  a  formal  approach, 
based  upon  techniques  of  estimation  theory.  Both  approaches  have  proven  useful  and,  in  many  cases,  the  ad  hoc 
synchronizers  fit  within  an  estimation  theory  framework.  We  follow  the  formal  approach,  deriving  a  structure  for 
joint  phase  and  timing  acquisition  and  tracking  based  upon  the  well-known  maximum  likelihood  (ML)  criterion  of 
optimality  (see  [4],  for  example)  which  maximizes  a  probability  density  function  for  received  samples  conditioned 
upon  a  set  of  parameters  to  be  estimated.  We  emphasize,  however,  that  because  of  necessary  approximations  for 
analytical  completion  and  receiver  implementation,  our  final  phase  and  timing  recovery  scheme  is  an  approximation  to 
the  ML  solution  and  is,  therefore,  suboptimal. 

Final  questions  involving  performance  of  this  suboptimal  synchronizer  will  require  simulation  methods  because  of 
inherent  analytical  complexities.  A  direct  measure  of  the  “goodness’"  of  an  estimator  is  its  variance,  usually  compared 
to  either  the  Cramer-Rao  or  the  Modified  Cramer-Rao  lower  bounds  (see  [  1  ]).  Ultimately,  desirable  performance  would 
be  in  terms  of  effects  on  ability  to  correctly  receive  data  when  using  the  particular  synchronization  method.  A  future 
paper  will  present  results  from  appropriate  simulation  studies. 
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II.  General  Structure  for  the  Phase  and  Timing  Synchronizer 

The  task  is  to  find  a  phase  and  timing  recovery  scheme  (acquisition  and  tracking)  that  can  operate  with  random  data 
and  known  frequency.  Furthermore,  the  synchronizer  is  to  be  non-data  aided ,  operating  without  benefit  of  a  known  data 
sequence  such  as  a  preamble,  and  non-decision  directed ,  operating  without  benefit  of  a  detected  data  sequence.  The 
notation  and  modeling  given  below  closely  follow  Mengali  and  D' Andrea  [1]  for  single-/;  CPM,  including  complex 
envelope  notation. 

We  assume  that  the  received  waveform,  r(t),  is  given  by 

r(t)  =  s(t)  +w{t),  (1) 

where 

=  T'Q)  (2) 

is  the  CPM  signal  and  wit)  is  an  additive  white  Gaussian  noise  component  with  spectral  density.  No-  In  equation  (2), 
0  and  r  represent  the  unknown  carrier  phase  and  time-of-arrival,  respectively;  Es  is  the  signal  energy  over  the  symbol 
interval,  T;  and  a  =  (. . . ,  a_i,  <*o,  «i,  •  •  ■)  is  the  data  sequence  with  <*,  €  (±1,  ±3}  for  the  4-ary  case  of  interest. 
Note  that  there  is  no  carrier  frequency  offset  since  carrier  frequency  is  assumed  to  be  known. 

The  phase,  ip(t,  a),  is  given  by 

-ip(t,a)  -  2Tr'^2aihiq{t  iT) ,  (3) 

i 

where  hi  is  the  modulation  index  for  the  ith  interval  and  q(t)  is  the  phase  response  function.  For  the  dual-/;  case, 
hi  =  ho  is  applied  to  c*o  =  (•••><*  2><To>Q;2>  •  •  •)  while  hi  =  hi  is  applied  to  ai  =  (. . . , a  3, a  1, <*1,0:3, . . .). 

Thus, 

ip(t,  a)  =  2nho  <*2 iq{t .....  2 iT)  +  2nh\  ^  <*2i+i<7(t  [2/  +  1]T)  (4) 

i  i 

-  ■fo(t,cxo)  +  fi(t,OCi)  . 

For  the  full-response,  rectangular  frequency  pulse  case  (so-called  1REC),  the  phase  response  is  a  ramp  function  given 
by 

{0,  if  t  <  0 

t/2T,  ifO  <t<T  (5) 

1/2,  if  t  >  T. 

In  keeping  with  modem  practice,  we  assume  a  digital  synchronizer  that  processes  discrete-time  samples  of  the 
waveforms  rather  than  an  analog  synchronizer  that  processes  directly  the  continuous-time  waveforms.  The  samples, 
x  =  {x(kTs}},  are  taken  as  shown  below  (Figure  1)  from  the  output  of  an  anti-aliasing  filter  (AAF),  which  is  assumed 
to  have  sufficiently  wide  passband  to  pass  the  signal,  s(t),  without  significant  distortion.  Also,  it  is  assumed  that  the 
sampling  period,  Ts,  is  sufficiently  short  to  allow  the  samples  to  retain  the  information  content  of  the  continuous-time 
waveforms. 


Fig.  1.  Sampling  output  of  anti-aliasing  filter  (AAF). 

Additionally,  we  let  Lq  be  the  number  of  information  symbols  to  be  “ observed ”  for  purposes  of  collecting  the 
samples  (Lq  is  assumed  to  be  an  even  integer),  and  let  N  be  the  number  of  samples  taken  per  symbol.  (Mengali  and 
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D ’Andrea  suggest  N  —  4  as  an  appropriate  value  for  cases  studied  by  them.)  Thus, 

To  =  “observation  time"  (6) 

=  L0T 
=  NL0TS 


and 

x  —  {x(0),x(Ts),  •  •  •  ,x([NLq  —  1]TS)}.  (7) 

The  task  is  to  derive  a  maximum-likelihood-based  estimator  of  the  unknown  phase  and  timing  which  does  not  require 
knowledge  of  a.  Hence,  the  likelihood  function,  A(x\9,t),  is  obtained  from  the  conditional  likelihood  function, 
A(:e|qi,  9,  f),  by  averaging  over  a  as  described  in  Mengali  and  D’ Andrea  [1].  Clarifying  previously  introduced  and 
additional  notation,  we  let 

9  —  the  actual,  but  unknown,  phase , 

9  —  a  possible  value  of  9  (a  “realization”), 

9  =  an  estimate  of  9, 

and  similarly  for  r,  r,  and  t. 

The  conditional  likelihood  function  is  given  by 


A(®|a,  9,  t)  =  exp  l  Re 


NL0  1 

Y  x(kTs)r(kTs ) 
L  fc=o 


^  §  “}■ 


where 


s(t)  =  ej\l‘^fei'i)(>(-t-T'OLo)ejxp i(t  T-Qi). 


Since  the  last  sum  in  equation  (8)  is  independent  of  a,  0.  and  r,  it  suffices  to  consider 


M{x\oc,9,t)  =  exp  l  jA  Re 


I  Ts  2  Es 

=  exp  fyVTRe 


iVL0  — 1 

Y  x(kTs)s*(kTs) 

.  k= 0 

NL0  1 


e  jd  Y  x(kT° 
k^O 


x  e  jv  o(kTs  T,(X0)e  jlpi(kTs  T,«i) 


(8) 


(9) 


(10) 


To  obtain  practical  results,  the  case  of  small  signal-to-noise  ratio  (SNR)  will  be  considered  since  a  synchronizer  that 
is  “optimum”  for  small  SNR  should  certainly  perform  well  with  higher  SNR.  For  small  SNR,  we  have 


A'(x\oc,9,t) 


=  A"(*|a,0,r), 


(ID 


where 


_  NLo  1 

Z  =  e~jd  Y  x(kTs)e  r,a0)e  H\{kTs  r,ai)_ 

k= 0 


Averaging  over  a  in  A"(®|a,  9,  r),  as  shown  in  Appendix  A,  yields 


(12) 


M\x\9,t) 


i+k\fWE«{Re[z]] 

1  +  !%)J  cos(0x  —  &)■> 


(13) 
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where 


and 


F(ho,hi,kTs  —  f) 


X  =  \X\e>*x  (14) 

NL0  1 

-  ^  x(kTs)F(ho,hi,kTs  r) 

A:=0 


JJcos  [ 2irhoq(kTs  t  2 iT)]  cos  [4 irhoq(kTs  r  2iT)]  (15) 

i 

x  cos  [2nhiq(kTs  —  f  —  [2i  +  1]T)]  cos  [4nhiq(kTs  —  r  —  [2 i  +  l]!1)] . 


Note  that  F(hoxhi,  kTs  —  t)  is  real- valued  and  independent  of  9. 

Next,  we  let  9 mi  and  9  ml  denote  the  maximum  likelihood  estimates  of  9  and  r,  respectively,  i.e.,  the  values  of  9  and 
r  that  maximize  A"(x\9,  9).  Since  |X|  is  independent  of  9.  it  follows  by  inspection  of  equation  (13)  that  9  ml  —  4>x> 
where  0X  is  defined  as  in  equation  (14).  (Technically,  9ml  and  9  ml  are  only  approximations  of  the  ML  estimates 
because  of  the  small  SNR  approximation  above.)  Note  that  9ml  depends  on  9  ml,  but  that  the  reverse  is  not  true. 
From  equations  (13)  and  (14),  it  follows  that  t ml  is  the  value  of  r  that  maximizes  |X|  or,  equivalently,  maximizes 
|X|2,  where 

NL0  INLo  1 

\X\2=  2  x(k1Ts)x*{k2Ts)F(h0,h1,klTs  T)F(ho,huk2Ts  r).  (16) 

ki=0  k2=0 

Since  it  is  mathematically  intractable  to  find  9 ml  as  a  solution  either  to  |X|  =  0  or  to  ^  |X|2  =  0,  we  will  describe 

a  practical  approach  to  determine  an  approximate  maximum  of  equation  (14)  and,  thus,  an  approximate  value  for  9  ml- 
We  divide  the  region  of  uncertainty  for  r,  namely,  the  time  interval,  T  —  NTS ,  into  Nr  increments,  each  of  width, 
At.  Letting  m  be  the  number  of  these  increments  contained  in  Ts,  it  follows  that 


Nt  = 


T_ 

Af 


NT8 

At 


Nm. 


(17) 


The  timing  synchronizer  may  be  implemented  as  shown  below  (Figure  2)  as  a  group  of  N  parallel  processors,  each  of 
which  generates  m  values  of  X  corresponding  to  m  assumed  values  of  f.  The  combination  produces  a  total  output  of 
Nm  values  of  X,  providing  the  basis  for  finding  9ml  and  9 ml  by  choosing  the  output  with  the  largest  modulus.  That 

is,  for  each  a  —  0, 1, ... ,  X in _ 1,  let  Ta  =  a  ■  At  be  the  ath  assumed  value  of  9  and  compute  the  corresponding 

value,  Xa,  by  substituting  Ta  into  equation  (14).  Next,  determine  the  value  oo  €  {0, 1, . . .  Nm  —  1}  such  that 


l-Xool  =  max{|Xa|  :  a  =  0, 1, . . . ,  JVm  -  1}. 


Finally,  it  follows  that 


9ml  =  Too  =  a0  •  At 


and,  since  Xao  =  \Xao  \  exp {j<l>x  ),  it  also  follows  that 


(18) 


0ml  -  4>xao  ■ 


(19) 
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Fig.  2.  Parallel  processors. 


Computation  of  Ajp  the  jth  output  of  the  ptJl  processor  (Figure  2),  for  p  —  1,2,. 

. . ,  N  and  j  =  0, 1, 2, . . . , 

m  __  1,  as  derived  in  Appendix  B,  is  performed  according  to  the  block  diagram  (Figure  3),  in 
and  Aj(fii),  are  defined  as 

which  the  scalars,  Kq(Jiq) 

Ko(ho)  =  cos(nho)  cos(2nho) 

(20) 

and 

N\{h\)  =  COs(7r/ll)  COs(27r/ll), 

(21) 

and  the  N  x  1  column  vectors,  Vq(/io,  m,j)  and  Vi(/ii,ra,  j),  are  given  by 


V0(h0,m,j) 


cos  [(7i7io/A0  (1  -  j/m)]  cos  [(27r/io/A^)  (1  —  j/m)] 
cos  [(7r/i0/Ar)  (2  -  j/m)]  cos  [(27iAo/A0  (2  -  j/m)] 

cos  [(7r/to/iV)  (TV  —  j/m)]  cos  [(2irho/N)  (AT  —  j/m)] 


and 


Vi(hi,m,j) 


cos  [(ithi/N)  (N  +  1  -  j/m)]  cos  [(27r/ii/iV)  (./V  +  1  —  j/m)] 
cos  [(nhi/N)  (N  +  2  -  j/m)]  cos  [(27r/z.1/iV)  (./V  +  2  —  j/m)] 

cos  [(7r/i,i/A^)  (2AT  —  j/m)]  cos  [(2nhi/N)  (2N  —  j/m)] 


The  lx  A  input  row  vectors,  xjj$,  are  defined  by 


=  [x{\p+(l-  1)N]  Ts),x(\p+(l-l)N  +  l]Ts),...,x(\p-l+  IN]  Ts)] , 


(22) 


(23) 


(24) 


f°r  l  —  1,2, ,  Lq;  for  those  values  of  indices  greater  than  A L{)  1,  the  upper  limit  for  the  actual  samples  given  in 

equation  (7),  we  set 


x(LoNTs)  =  x([L0N  +  1]TS)  =  •  •  •  =  x{[LqN  +  p-  1  ]TS)  =  0. 


(25) 
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For  each  fixed  value  of  p,  m  different  values,  { X ^  i)m+j  :  j  =  0, 1, 2, . . . ,  m  —  1},  are  generated  corresponding  to 
the  m  different  values  of  the  vector  pair,  Vo(ho,  m,j)  and  Vi(/ii,ra,  j),  Asp  varies,  the  input  row  vectors  change  as 
does  the  input  to  the  final  summation  producing  X(p_i)m+J-. 


*:  x 


H  Delay  2N  Samples 


X 

* — ( x — 


—  x  } — 
T 

V0(hQ.mJ) 


Jp)  Jp) 

\Lo-W . X1A'  / 


Output  of  Group  p 
Organizer 


<A> / 2)— I 
1=0 


x  [£,,(*>)*, 


ATC  seconds 


*([p  i]7; x(ry).jr(0) 


x. 


V2 

A'o(*o>Z  xVNv\{h\-m-J) 

v(  P)  v(P)  A 

*L0NT"'2N 

V /-I 

x  [AWAA^r1 

X 


*o(W(*i) 


Delay  2Ar  Samples 


f 


Fig.  3.  Calculation  of  X(p_i)m+J  by  the  pth  parallel  processor,  for  p  =  1, 2, . . . ,  Ar  and  j  =  0, 1, . . . ,  m  1. 

Summarizing,  the  calculation  of  X(p  and  thus  off  ml  and  0  ml*  is  performed  as  described  in  the  following 

algorithm. 

Step  1:  Pass  the  received  signal,  r(t),  through  the  AAF,  and  sample  the  resulting  band-limited  waveform  (Figure  1). 

Step  2:  Collect  and  organize  the  NLq  samples  into  N  groups  (each  group  includes  all  NLq  samples)  for  input  into 
N  sets  of  parallel  processors  as  described  below. 

Group  1:  x(0),x$,  x^, . . .  Xl0n- 

Group  2:  x(0),  x(Ts), x($,  x$, . . . ,  xJJ^. 

Group p:  x(0),  x(Ta),  x(2Ts), . . . ,  x(\p  1]TS),  x£$,  x^, . . . ,  x^. 

Group N:  x{0) , x(Ts) , x(2Ts) , . . . ,  x([N  1  }TS), x^\ x^}, . . . , x^. 

Step  3:  Calculate  the  Nm  values  of  X,  using  N  parallel  sets  of  processors.  For  example,  the  pth  set  of  processors 
associated  with  Group  p  input  above,  produces  the  m  values,  {  A^,  \)m+j  '■  3  —  0, 1, 2, . . . ,  m  —  1}. 

Step  4:  Compute  {\X(p  i  )m+j  |2  :  p  —  1, 2  j  =  0, 1, 2, ...  ,m  —  1}  and  choose  the  largest  value,  say 

\Xa,} |2,  where  Xao  =  \Xao\exp(j<j>XaJ. 

Step  5:  Compute  tml  and  0ml ,  as  in  equations  (18)  and  (19),  respectively. 


Preceding  Page  Blank 
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III.  Compact  Description  of  Phase  and  Timing  Synchronizer 

From  the  NLq  samples,  x(0),x(Ts), ...,  x([NLo  l]Ta),  form  the  row  vectors  x\^j  previously  defined  in  equa¬ 
tion  (24)  for p  —  1,2, ...  ,N  and  l  —  1,2,..., Lq.  We  note  from  equations  (24)  and  (25),  when  l  =  Lq  and  p  =  N, 
then 

XS  =  [°>  0,...,0],  (26) 

and  we  extend  the  definition  of  x|^  to  the  case  when  1  =  0  and  p  =  1  by  defining  the  vector 

4n  =  <T.), . . . ,  x([N  1]TS)] .  (27) 

We  also  define  the  1  x  N  row  vector 

P(p)  =  [1,1,..  .,1,0,0,  ...0J  (28) 

P  A'  p 

and  the  1  x  Lq/2  row  vector 


K0  = 


1,  Ko(ho)Ki(hi),  {Ko^KM}2  ,...,  {A'o(M^iOi)} 


1 

2  1 


(29) 


where  KqUiq),  Ki(h\)  are  as  defined  previously  in  equations  (20)  and  (21),  respectively.  Then,  we  define  the  N  x  m 
matrices 

V(°)  =  [V0(/io,  m ,  0),  Vo(h0,  m,  1), . . . ,  V0(/io ,  m,m-  1)]  (30) 

and 

v(1)  =  [Vi(hi,m,0),V1(hl,m,l),. . .  -  1)],  (31) 

where  Vo(ho,m,  j)  and  Ni(h\ ,m,j)  were  defined  previously  in  equations  (22)  and  (23),  respectively.  Next,  letting 
AT  denote  the  transpose  of  A,  we  define  three  other  matrices  that  are  each  functions  of  the  NLq  samples:  the  N  x  m 
matrix. 


XOA 


(4a)  >(4a 


m 


(32) 


and  the  Lq /2  x  N  matrices, 


,(p)  __ 
odd 


and 


x(p)  = 


Jp) 
'IN 
C(P) 
c3  N 


Y(P) 

L  X(L0  1)A 


((P) 

t2Ar 

,(p) 

Hn 


(33) 


„(p) 

L  XI0A  J 


(34) 


Finally,  using  the  vectors  and  matrices  just  defined,  the  calculation  of  the  Nm  values  of  X,  namely  Xq,  X\, ... , 
Xjvm  i,  can  be  represented  in  compact  form  in  terms  of  matrix  operations.  As  shown  in  the  block  diagram  (Figure  3), 
a  general  value  of  X,  say  X(p  i )m+j,  is  computed  as  the  grand  total  of  two  summation  expressions  in  the  upper  and 
lower  portion  of  the  diagram  and  a  sum  of  sample  values  in  the  middle.  For  a  fixed  value  of  p  and  with  j  varying 
between  0  and  m  _  1,  it  can  be  readily  shown  that  Kox^V^  ant*  Ao  ( /i  o )  K  o  xlP  e  n  V  *■ 1  represent  the  values  of  the 
summation  expressions  in  the  upper  and  lower  portions  of  the  block  diagram,  respectively,  and  that  P^xoa  represents 
the  sum  of  sample  values  in  the  middle.  Hence,  all  Nm  values  of  X  can  be  calculated  as 

\p^-(p  J)mi  •  •  •  ;  •  •  •  >  Apm  \  ]  =  P^XOA  +  Kox^V(°)  +  AotMKox^vW,  (35) 

forp  =  1,2,...,  A. 
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IV.  Conclusions 

This  document  presents  the  derivation  of  a  joint  phase  and  timing  synchronizer,  assuming  frequency  is  known, 
for  dual-//,  full-response,  4-ary  CPM  with  1REC  instantaneous  frequency  pulse.  The  synchronizer  is  based  upon  the 
maximum-likelihood  criterion,  with  the  final  structure  resulting  from  certain  approximations  made  to  overcome  various 
mathematical  complexities  and  to  make  any  subsequent  implementation  of  reasonable  complexity. 

It  is  not  clear  that  this  new  synchronizer  provides  significant,  if  any,  benefits  compared  with  others  in  use  or  proposed 
for  specific  applications.  Comparisons  based  upon  simulation  are  required  before  any  further  conclusions  can  be  made. 
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Appendix  A 

Conditional  Likelihood  Function  Averaged  Over  Random  Data  Symbols 

The  conditional  likelihood,  A"(x\6,  f),  given  in  equation  (13),  is  obtained  by  averaging  R e[Z]  over  a,  with  Z  given 
in  equation  (12).  Since  the  data  symbols  occur  independently,  we  consider  separately  the  two  averages, 

Ea0  (exP  [~j%(kTs  ~  t.  a0)]} 


E, 

Sl  {exp {-jipy/kTs  -r.a])]}. 

Consider  the  first  term: 

£a-.  W’  Hv0(*T.  -  t,50)]  1  =  | 

|exp 

j2nho^ra2iq(kTs  r  2  iT) 

i 

=  £a.| 

[ilexP  [  -j2irhooi2iq{kT8  t  2TT)]j 

=  Ea2t  {exp  [-j2nhoa2iq{kTs  -  f-  2 iT)}} 

l 


—  JJ(l/4)  {exp  [j&nhoq(kTs  —  f  —  2 iT)]  +  exp  {—j§jnhoq{kTs  —  r  —  2 iT)) 

i 

+  exp  [j2'nhaq{kTs  —  f  —  2  iT)]  +  exp  [—j2'Khoq(kTs  —  r  —  2  iT)]} 

—  JJ(l/2)  {cos  [Qnhoq(kTs  -  t  -  2 iT)]  T  cos  [2nhoq(kTs  —  r  —  2 iT)]} . 

i 

Similarly, 

E0Ci  {exP[-j'V'i(fcT5  -  T,ai)]}  =  J"J(l/2)  {cos  [Qnhlq(kTs  -  t  -  [2i  +  1}T)] 

i 

+  cos  [2nhiq(kTs  —  r  —  [2i  +  l]r)]}  . 

Using  the  identity,  cos  3/3  +  cos  (3  —  2  cos  3  cos  2/3,  yields  the  result  in  equation  (13). 
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Appendix  B 

Calculation  of  X^p  1)m+j 

In  this  appendix  we  derive  the  mathematical  expression  for  X(p  \)+j*  the  output  of  the  parallel  processor  in  the 
general  case  (Figure  3).  Based  on  the  definition  of  ra  (section  II)  and  equation  (14),  it  follows  that 


NL0  1 

X{p . l)m+j  —  E  x(kTs)F(ho,hi,kTs  ra), 

k= o 

where  fa  =  a-  At  and  a  =  (p  —  1  )m+j,  p—  1,2,...,  N,  and  j  =  0, 1, . . . ,  m  —  1.  Expanding  the  above  summation 
into  several  sums  yields 


(p  1  )m+j  — 


pi  (p  1)+Ar 

^2x(kTs)F(h0,hi,  kTs  -  fa)  +  ^  x(kT$)F(h0,  hi,  kTs  -  fa)  +  . . . 

k= 0  k—p 

(p.A)+rN  (p  l)+L0N 

...+  x(kTs)F(ho,hi,k,Ts  -  fa)  4-  . . .  4-  ^  x(kTs)F(ho,  hi,  kTs  -  ra) 

k=p+(r  .  1)N  k=p+{L0  l)N 

pl  L0  (P  l)+r Ar 

=  ^x(fcTs)F(/i0,hi,A;rs  ~Ta)  4-^  ^  x(kTs)F{h0,hi,kTs -Ta),  (B.l) 

fc=0  r=  1  k=p+(r  V)N 


where  F (ho,  hi,  kTs  —  f  a)  is  the  product  defined  in  equation  ( 1 5)  with  indices  i  =  0, 1, . . . ,  Lq/2  —  1.  (We  remark  that 
the  last  p  —  1  values  of  x(kTs)  in  the  second  summation  in  the  last  line  of  equation  (B.l)  are  0,  as  noted  in  equation  (25) 
above.) 

The  function  F(ho,hi,kTs  _  r)  is  very  complicated,  having  Lq/2  sets  of  four  cosine  factors  and  exhibiting  a 
complex  interaction  between  the  summation  index  k  of  equation  (14)  and  the  product  index  i  of  equation  (15).  The 
calculation  of  each  X(p  xnn+j  requires  repeated  calculation  of  F(ho,  h\ ,  kTs  fa),  which  can  be  broken  down  into 
four  cases,  depending  on  the  values  of  fa,  k,  and  i  as  described  below. 

Case  1:  (k  _  2 Ni)Ts  <  fa 

cos[27r  hoq(kTs  —  Ta  —  2  iT)\  =  cos  [Anhoq{kTs  —  fa  —  2  iT)\  =  1 
cos  [2iThiq(kTs  —  ra  —  [2  i  4-  1]T)]  =  cos  [4irhiq(kTs  —  fa  —  [2  i  +  1]T)]  =  1 

Case  2:  ra  <  {k  —  2 Ni)Ts  <T  +  ra 


cos  [2nhoq(kTs  —  fa  —  2iT)\ 
cos  [Anhoq(kTs  —  fa  —  2  iT)\ 
cos  [2nhiq(kTs  —  ra  —  [2  i  +  1]T)] 

Case  3:  T  4-  fa  <  (k  —  2 Ni)Ts  <2 T  +  rQ 

cos[2Trhoq(kTs  —  ra  —  2iT)\ 
cos  [Anhoq{kTs  —fa  —  2  iT)\ 
cos  [2nhiq(kTs  —  fa  —  [2  i  +  1]T)] 
cos  [4irhiq(kTs  —  ra  —  [2  i  4-  1]T)] 

Case  4:  ( k  —  2Ni)Ts  >  2T  +  fa 


cos  [27rho  {(k/2N)  -  i-  ( fa/2NTs )}] 
cos  [47r/io  {( k/2N )  -  i  -  ( fa/2NTs }] 
cos  [Airhiq(kTs  —  fa  —  [2 i  4-  1]T)]  =  1 

COS  7r/lo 

cos  2irho  ® 

cos  [27rhi  {(k/2N)  -  i  -  ( fa/2NTs )}] 

cos  [4tt hx  {(k/2N)  -  i  -  (7a/2NTs)}}  ^ 

c 

o 


cos[27rhoq(A:Ts  — 
cos  [Anhoq(kTs  —fa  —  2iT)\  =  cos  2tt/io 
cos  [2nhiq(kTs  —  fa  —  [2i  4-  1]T)]  =  cos  n  hi 
cos  [4nhiq(kTs  —  fa  —  [2z  4-  1]T)]  =  cos  2nhi 
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Determining  which  of  the  four  cases  applies  for  a  particular  (i,  k)  pair  is  greatly  facilitated  by  drawing  in  the  (i,  k)- 
plane  (Figure  4)  the  boundaries  k  =  2 Ni,  k  =  (2i  +  1  )N,  and  k  =  (2i  +  2)N,  which  apply  for  r  =  0,  and  then 
inferring  the  effects  of  f  ^  0  on  the  choice  of  factors. 


k 


Fig.  4.  Regions  for  factors  of  F(ho,  hi ,  kTs  r). 

For  k  =  0, 1, . . .  ,p  1  and  i  —  0, 1, . . . ,  Lo/2  1,  it  follows  that 

(k  _  2Ni)Ts  —  ( k  2Ni)m  ■  Af 

<  [(p  —  1)  —  2Ni ]  m  ■  Af 

—  [( p  —  1  )m  —  2 Nim]  ■  Af 

^  [ip  ~  +  j)  ■  Ar 

—  a  ■  Af 
=  Ta, 

where  the  first  equality  holds  by  the  definition  of  Ts  in  equation  ( 1 7),  the  first  inequality  holds  since  0  <  k  <  p  —  1, 
and  the  second  inequality  holds  since  2Nim  >  0  and  j  >  0.  Thus,  Case  1  holds  for  all  values  of  A:  in  the  first  sum  in 
equation  (B.l)  above,  which  implies  that  F(ho,  h\,kTs  -  fa)  =  1  for  such  k  and  that,  therefore,  this  sum  is  equal  to 
X^=o  x(kTs).  Hence,  it  follows  that 


P— 1  Lo  (p~l)+rN 

A(p  1  )m+j  =  E  x(kTs)  +  E  E  x(kTs)F(ho,hi,kTs  fa ), 

k= 0  r= 1  k-p+(r  1)N 

which  can  be  expanded  as 

P  1  L0  ip  l)+rN  L0  (P  1)+WV 

X(jp  1  )m+j  —  J2  x(kTs)+J2  E  *(KTe)F(ho,  hu  kT8  ra)+  E  E  x(kTs)F(h0,  hi, kTs— Ta). 

k=0  r=l  k=p+(r._l)N  r=l  k=p+(r  1)N 

r  odd  r  even 

Noting  that  the  odd  values  of  r  between  1  and  Lo  may  be  expressed  as  r  —  21  +  1,  for  l  =  0, 1, . . . ,  Lo/2  _  1,  and  the 
even  values  of  r  may  be  expressed  as  r  =  21,  for  l  -  1, 2, . . . ,  L0/2,  we  observe  that  the  last  equation  for  X(p_i)m+j 
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may  be  written  equivalently  as 


p- 1  (Lo/2)  1  (p.  l)+(2l+l)N 

E  x(kTs)  +  E  E  x(kTs)F(h0,hi,kTs  ra) 

k—0  1=0  k=p+2lN 

L0/2  (p.  1)+2W 

+  E  E  x(kTs)F(h0,hi,kTs  ra).  (B.2) 

i=l  k=p+(2l  1)N 

To  verify  the  calculation  of  Xip  j  )m .  j  depicted  in  the  block  diagram  above  (Figure  3),  we  observe  that  the  value  of 
the  first  sum  in  equation  (B.2)  equals  the  sum  of  sample  values,  x(\p  l]Tj), . . . ,  x(Ts),  x(0),  entering  the  summation 
operator  in  the  middle  of  the  figure.  We  claim  further  that  the  values  of  the  second  and  third  sums  in  equation  (B.2) 
equal  the  values  output  by  the  upper  and  lower  blocks,  respectively,  which  enter  the  summation  operation  in  the  figure. 
We  now  prove  these  claims. 

To  evaluate  the  two  double  sums  in  equation  (B.2),  the  interaction  between  the  summation  index,  k,  of  the  inner 
sums  and  the  product  index,  i,  of  equation  (15)  must  be  understood.  After  some  analysis,  the  combinations  of  k  and 
i  belonging  to  each  of  the  four  cases  described  above  can  be  determined  as  tabulated  (Table  I).  Using  the  results  for 
k  =  p  +  2 IN  +  b  to  evaluate  F(ho,  h\,  kTs  fa),  substituting  [ip  1  )rn  +  j]Ts/m  for  ra,  and  making  an  obvious 
change  of  summation  index  yields  that 

(p  l)+(2l+l)N  yv  1 

^  x(kTs)F(ho,  hi,  kTs  -  Ta)  =  [/sT0(^o)^i  {hi)]1  Y  x  (l-P  +  2lN  +  b] T*) 

k=p+2lN  6=0 

x  cos  \(irho/N){\  +  b  —  j/m )] 
x  cos  [(27rho/A0(l  +  b  —  j/m)} 

=  [Ko(ho)K\  (hi)]/  x^+1)ArV0(ho,  m,j),  (B.3) 

where  Ko(ho ),  K\  (hi),  x[^+1)/v  and  Vo(ho,m,j)  are  as  defined  in  equations  (20),  (21),  (24),  and  (22),  respectively. 
Similarly,  using  the  results  for  k  =  p  +  (21  l)N  +  b.  it  can  be  shown  that 

(p  l)+2 IN  N  1 

x(kTs)F(ho,hx,kTs -To)  =  Ko^iKoihomihi)}1  1  ^x([p+(2(  A)N  +  b}Ts) 

k=p+(2l  1  )N  6=0 

x  cos  [(nhi/N)(N  +  1  4-  b  —  j/m)] 
x  cos  [(27rhi/Ar)(iV  +  1  +  h  —  j/m)] 

=  Ko(h0)  [Ko(h0)Ki  (hi)f  1  x^- V i (h\ ,m,j),  (B.4) 

where  and  V  i  (hi,  m,j)  are  as  defined  in  equations  (24)  and  (23),  respectively.  Substitution  of  equations  (B.3) 
and  (B.4)  into  equation  (B.2)  proves  the  two  claims  above,  which  completes  the  derivation  of 

TABLE  I 

Values  of  i  satisfying  cases  L  2.  3,  and  4,  for  given  values  of  k  (b  =  0, 1, . . . ,  N  1). 


Values  of  i  Satisfying: 

1 

k 

Case  1 

Case  2 

Case  3 

Case  4 

0,  l,...,(Lo/2)_l 

l,2,...,L0/2 

P  +  2lN  +  b 
p+(2l-l)N  +  b 

(*f  1,1 +2,...,  (Lo/2)  1 

Lq/2 

l 

None 

None 

FI 

0,1,...,Z_1 
0,1,...,/ _2 
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